%%%�����Է���
%% date:2018-08-24
%% by Yaokui Cui 
% From Institute of Remote Sensing and GIS,School of Earth and Space Sciences, Peking University
% Email:yaokuicui@pku.edu.cn
% 2021-02-18
% delete(gcp)
% parpool
clear
years={'2002';'2003';'2004';'2005';'2006';'2007';'2008';'2009';'2010';'2011';'2012';'2013';'2014';'2015'};
for i_str=12:12
    for i_year=2012:2012
        mkdir(['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_e\' ,num2str(i_year)]);
        mkdir(['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_t\' ,num2str(i_year)]);
        list_albedo=dir(['G:\CAS\data\ICTP\ET_Heihe\albedo\',num2str(i_year),'\','*.tif']);
        list_e=dir(['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\E_RS_heihe\',num2str(i_year),'\','*.tif']);
        list_t=dir(['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\T_RS_heihe\',num2str(i_year),'\','*.tif']);
        list_lrad=dir(['G:\CAS\data\ICTP\ET_Heihe\LRad\',num2str(i_year),'\','*.tif']);
        list_ndvi=dir(['G:\CAS\data\ICTP\ET_Heihe\NDVI\',num2str(i_year),'\','*.tif']);
        list_pres=dir(['G:\CAS\data\ICTP\ET_Heihe\Pres\',num2str(i_year),'\','*.tif']);
        list_sm1=dir(['G:\CAS\data\ICTP\ET_Heihe\SoilMoisture_ASCAT2_1km5_noLL\1km_daily\',num2str(i_year),'\','*.tif']);
        list_srad=dir(['G:\CAS\data\ICTP\ET_Heihe\SRad\',num2str(i_year),'\','*.tif']);
        list_temp=dir(['G:\CAS\data\ICTP\ET_Heihe\Temp\',num2str(i_year),'\','*.tif']);
        list_rh=dir(['G:\CAS\data\ICTP\ET_Heihe\rh\',num2str(i_year),'\','*.tif']);
        list_wind=dir(['G:\CAS\data\ICTP\ET_Heihe\wind\',num2str(i_year),'\','*.tif']);
        DEM=imread('G:\CAS\data\ICTP\ET_Heihe\DEM.tif');
        DEM=single(DEM);
        DEM=(DEM<0).*0+(DEM>0).*DEM;
        landcover=imread('G:\CAS\data\ICTP\ET_Heihe\Landcover2009.tif');
        k=365;
        R_result=zeros(550,600);
        %   lst_365=zeros(550,600,365);
        ndvi_365=zeros(550,600,365);
        %     lai_365=zeros(550,600,365);
        albedo_365=zeros(550,600,365);
        wind_365=zeros(550,600,365);
        temp_365=zeros(550,600,365);
        rh_365=zeros(550,600,365);
        pres_365=zeros(550,600,365);
        lrad_365=zeros(550,600,365);
        srad_365=zeros(550,600,365);
        e_365=zeros(550,600,365);
        t_365=zeros(550,600,365);
        lc_365=zeros(550,600,365);
        dem_365=zeros(550,600,365);
        ET_25km_result_out=uint8(zeros(550,600,365));
        %%%%import data
        for i=1:365
            dem_365(:,:,i+(i_year-2012)*365)=DEM;
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\NDVI\',num2str(i_year),'\', list_ndvi(i).name);
            ndvi_25km=imread(str);
            ndvi_365(:,:,i+(i_year-2012)*365)=ndvi_25km;
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\albedo\',num2str(i_year),'\', list_albedo(i).name);
            albedo_25km=imread(str);
            albedo_365(:,:,i+(i_year-2012)*365)=albedo_25km;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\wind\',num2str(i_year),'\', list_wind(i).name);
            temp=imread(str);
            wind_365(:,:,i+(i_year-2012)*365)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\temp\',num2str(i_year),'\', list_temp(i).name);
            temp=imread(str);
            temp_365(:,:,i+(i_year-2012)*365)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\rh\',num2str(i_year),'\', list_rh(i).name);
            temp=imread(str);
            rh_365(:,:,i+(i_year-2012)*365)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\pres\',num2str(i_year),'\', list_pres(i).name);
            temp=imread(str);
            pres_365(:,:,i+(i_year-2012)*365)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\lrad\',num2str(i_year),'\', list_lrad(i).name);
            temp=imread(str);
            lrad_365(:,:,i+(i_year-2012)*365)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\srad\',num2str(i_year),'\', list_srad(i).name);
            temp=imread(str);
            srad_365(:,:,i+(i_year-2012)*365)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\SoilMoisture_ASCAT2_1km5_noLL\1km_daily\',num2str(i_year),'\', list_sm1(i).name);
            temp=imread(str);
            sm1_365(:,:,i+(i_year-2012)*365)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\E_RS_heihe\',num2str(i_year),'\', list_e(i).name);
            temp=imread(str);
            e_365(:,:,i+(i_year-2012)*365)=temp;
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\T_RS_heihe\',num2str(i_year),'\', list_t(i).name);
            temp=imread(str);
            t_365(:,:,i+(i_year-2012)*365)=temp;
            lc_365(:,:,i+(i_year-2012)*365)=landcover;
        end
    end
    
    albedo_365=albedo_365/10000.;
    albedo_365=(albedo_365>0).*(albedo_365<0.3).*albedo_365;
    ndvi_365=ndvi_365/10000.;
    ndvi_365=(ndvi_365>0.1).*(ndvi_365<=1).*ndvi_365;
    
    for j_ns=1:1%%6:12%4:8%1:15 ;;int((55-40))/2+1:int((55-25))/2+1+1:;;ceil jin yi fa
        for j_nl=1:1 %%1:13%1:40����int((70-70))/2+1:int((105-70))/2+1+1:
            m=1;
            SM_25km_result=zeros(550,600);
            albedo_all=zeros(1,50000);
            ndvi_all=zeros(1,50000);
            wind_all=zeros(1,50000);
            temp_all=zeros(1,50000);
            rh_all=zeros(1,50000);
            pres_all=zeros(1,50000);
            srad_all=zeros(1,50000);
            lrad_all=zeros(1,50000);
            sm1_all=zeros(1,50000);
            dem_all=zeros(1,50000);
            e_all=zeros(1,50000);
            t_all=zeros(1,50000);
            lc_all=zeros(1,50000);
            
            for i=1:365
                ndvi_25km=ndvi_365(:,:,i);
                albedo_25km=albedo_365(:,:,i);
                wind_25km=wind_365(:,:,i);
                temp_25km=temp_365(:,:,i);
                rh_25km=rh_365(:,:,i);
                pres_25km=pres_365(:,:,i);
                srad_25km=srad_365(:,:,i);
                lrad_25km=lrad_365(:,:,i);
                sm1_25km=sm1_365(:,:,i);
                dem_25km=dem_365(:,:,i);
                e_25km=e_365(:,:,i);
                t_25km=t_365(:,:,i);
                lc_25km=lc_365(:,:,i);
                c=e_25km>0.1&t_25km>0.1&e_25km<7&t_25km<7&ndvi_25km>0.011&albedo_25km>0.011 &lc_25km >0 &lc_25km<16;
                c1=numel(ndvi_25km(c));
                
                if c1==0
                    continue
                end
                if c1>0
                    albedo_all(m:m+c1-1)=albedo_25km(c);
                    ndvi_all(m:m+c1-1)=ndvi_25km(c);
                    wind_all(m:m+c1-1)=wind_25km(c);
                    temp_all(m:m+c1-1)=temp_25km(c);
                    rh_all(m:m+c1-1)=rh_25km(c);
                    pres_all(m:m+c1-1)=pres_25km(c);
                    srad_all(m:m+c1-1)=srad_25km(c);
                    lrad_all(m:m+c1-1)=lrad_25km(c);
                    sm1_all(m:m+c1-1)=sm1_25km(c);
                    dem_all(m:m+c1-1)=dem_25km(c);
                    e_all(m:m+c1-1)=e_25km(c);
                    t_all(m:m+c1-1)=t_25km(c);
                    lc_all(m:m+c1-1)=lc_25km(c);
                    m=m+c1;
                end
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%train%%%%%%%%%%%%
    m
    c=e_all>0&t_all>0;
    
    albedo_all_train=albedo_all(c);
    [albedo_all_train,PALBEDO]=mapminmax(albedo_all_train,0.011,0.95);
    
    ndvi_all_train=ndvi_all(c);
    [ndvi_all_train,PNDVI]=mapminmax(ndvi_all_train,0.011,0.95);
    
    wind_all_train=wind_all(c);
    [wind_all_train,Pwind]=mapminmax(wind_all_train,0.011,0.95);
    
    temp_all_train=temp_all(c);
    [temp_all_train,Ptemp]=mapminmax(temp_all_train,0.011,0.95);
    
    rh_all_train=rh_all(c);
    [rh_all_train,Prh]=mapminmax(rh_all_train,0.011,0.95);
    
    pres_all_train=pres_all(c);
    [pres_all_train,Ppres]=mapminmax(pres_all_train,0.011,0.95);
    
    srad_all_train=srad_all(c);
    [srad_all_train,Psrad]=mapminmax(srad_all_train,0.011,0.95);
    
    lrad_all_train=lrad_all(c);
    [lrad_all_train,Plrad]=mapminmax(lrad_all_train,0.011,0.95);
    
    sm1_all_train=sm1_all(c);
    [sm1_all_train,Psm1]=mapminmax(sm1_all_train,0.011,0.95);
    
    
    lc_all_train=lc_all(c);
    [lc_all_train,Plc]=mapminmax(lc_all_train,0.011,0.95);
    
    dem_all_train=dem_all(c);
    [dem_all_train,Pdem]=mapminmax(dem_all_train,0.011,0.95);
    
    P=[ndvi_all_train(:) albedo_all_train(:) wind_all_train(:) temp_all_train(:) rh_all_train(:) ...
        pres_all_train(:) srad_all_train(:) lrad_all_train(:) dem_all_train(:) sm1_all_train(:) lc_all_train(:) ....
        ];
    P=P';
    e_all_train=e_all(c);
    [e_all_train,Pe]=mapminmax(e_all_train,0.011,0.95);
    
    t_all_train=t_all(c);
    [t_all_train,Pt]=mapminmax(t_all_train,0.011,0.95);
    
    
    T=[e_all_train(:) t_all_train(:)];
    T=T';
    
    inputSize = 11;
    
    layers = [ ...
        sequenceInputLayer(inputSize) %�����
        fullyConnectedLayer(128)
        reluLayer
        % eluLayer
        fullyConnectedLayer(128)
        reluLayer
        
        %     tanhLayer
        %      fullyConnectedLayer(16)
        %     reluLayer
        %     fullyConnectedLayer(8)
        %     reluLayer
        fullyConnectedLayer(2)
        %     softmaxLayer %softmax����
        % dropoutLayer(0.2)
        regressionLayer
        ]
    
    idx = randperm(size(P,2),300000);
    PValidation = P(:,idx);
    P(:,idx) = [];
    TValidation = T(:,idx);
    T(:,idx) = [];
    
    options = trainingOptions('adam', ...
        'ValidationData',{PValidation,TValidation}, ...
        'SequenceLength','shortest',...
        'ValidationFrequency',10, ...
        'LearnRateSchedule','piecewise', ...
        'InitialLearnRate',0.05, ...
        'LearnRateDropFactor',0.9, ...
        'LearnRateDropPeriod',200, ...
        'ExecutionEnvironment','gpu', ...
        'GradientThreshold',1, ...
        'MaxEpochs',1000, ...
        'MiniBatchSize',512, ...
        'Plots','training-progress');
    
    net = trainNetwork(P,T,layers,options);
    SM_all_out = predict(net,P);
    
    a=corrcoef(T,SM_all_out);
    a(1,2)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
    for i_year=2012:2012
        mkdir(['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_e\' ,num2str(i_year)]);
        mkdir(['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_t\' ,num2str(i_year)]);
        mkdir(['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_et\' ,num2str(i_year)]);
        list_albedo=dir(['G:\CAS\data\ICTP\ET_Heihe\albedo\',num2str(i_year),'\','*.tif']);
        list_e=dir(['G:\CAS\data\ICTP\ET_DNN_TSEB\E_RS\',num2str(i_year),'\','*.tif']);
        list_lrad=dir(['G:\CAS\data\ICTP\ET_Heihe\LRad\',num2str(i_year),'\','*.tif']);
        list_ndvi=dir(['G:\CAS\data\ICTP\ET_Heihe\NDVI\',num2str(i_year),'\','*.tif']);
        list_pres=dir(['G:\CAS\data\ICTP\ET_Heihe\Pres\',num2str(i_year),'\','*.tif']);
        list_sm1=dir(['G:\CAS\data\ICTP\ET_Heihe\SoilMoisture_ASCAT2_1km5_noLL\1km_daily\',num2str(i_year),'\','*.tif']);
        list_srad=dir(['G:\CAS\data\ICTP\ET_Heihe\SRad\',num2str(i_year),'\','*.tif']);
        list_temp=dir(['G:\CAS\data\ICTP\ET_Heihe\Temp\',num2str(i_year),'\','*.tif']);
        list_rh=dir(['G:\CAS\data\ICTP\ET_Heihe\rh\',num2str(i_year),'\','*.tif']);
        list_wind=dir(['G:\CAS\data\ICTP\ET_Heihe\wind\',num2str(i_year),'\','*.tif']);
        DEM=imread('G:\CAS\data\ICTP\ET_Heihe\DEM.tif');
        DEM=single(DEM);
        DEM=(DEM<0).*0+(DEM>0).*DEM;
        landcover=imread('G:\CAS\data\ICTP\ET_Heihe\Landcover2009.tif');
        k=365;
        R_result=zeros(550,600);
        ndvi_365=zeros(550,600,365);
        lai_365=zeros(550,600,365);
        albedo_365=zeros(550,600,365);
        wind_365=zeros(550,600,365);
        temp_365=zeros(550,600,365);
        rh_365=zeros(550,600,365);
        pres_365=zeros(550,600,365);
        lrad_365=zeros(550,600,365);
        srad_365=zeros(550,600,365);
        rn_365=zeros(550,600,365);
        sm1_365=zeros(550,600,365);
        sm2_365=zeros(550,600,365);
        e_365=zeros(550,600,365);
        t_365=zeros(550,600,365);
        lc_365=zeros(550,600,365);
        dem_365=zeros(550,600,365);
        ET_25km_result_out=uint8(zeros(550,600,365));
        %%%%import data
        for i=1:365
            dem_365(:,:,i)=DEM;
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\NDVI\',num2str(i_year),'\', list_ndvi(i).name);
            ndvi_25km=imread(str);
            ndvi_365(:,:,i)=ndvi_25km;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\albedo\',num2str(i_year),'\', list_albedo(i).name);
            albedo_25km=imread(str);
            albedo_365(:,:,i)=albedo_25km;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\wind\',num2str(i_year),'\', list_wind(i).name);
            temp=imread(str);
            wind_365(:,:,i)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\temp\',num2str(i_year),'\', list_temp(i).name);
            temp=imread(str);
            temp_365(:,:,i)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\rh\',num2str(i_year),'\', list_rh(i).name);
            temp=imread(str);
            rh_365(:,:,i)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\pres\',num2str(i_year),'\', list_pres(i).name);
            temp=imread(str);
            pres_365(:,:,i)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\lrad\',num2str(i_year),'\', list_lrad(i).name);
            temp=imread(str);
            lrad_365(:,:,i)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\srad\',num2str(i_year),'\', list_srad(i).name);
            temp=imread(str);
            srad_365(:,:,i)=temp;
            
            str= strcat ('G:\CAS\data\ICTP\ET_Heihe\SoilMoisture_ASCAT2_1km5_noLL\1km_daily\',num2str(i_year),'\', list_sm1(i).name);
            temp=imread(str);
            sm1_365(:,:,i)=temp;
            
            lc_365(:,:,i)=landcover;
            
        end
        albedo_365=albedo_365/10000.;
        lai_365=lai_365*0.1;
        lai_365=(lai_365>0.1).*(lai_365<10).*lai_365;
        
        albedo_365=(albedo_365>0).*(albedo_365<0.3).*albedo_365;
        ndvi_365=ndvi_365/10000.;
        ndvi_365=(ndvi_365>0.1).*(ndvi_365<=1).*ndvi_365;
        for i=1:365
            
            SM_25km_result=zeros(550,600);
            SM_25km_e=zeros(550,600);
            SM_25km_t=zeros(550,600);
            SM_25km_et=zeros(550,600);
            %                     lst_25km=lst_365(:,:,i);
            ndvi_25km=ndvi_365(:,:,i);
            lai_25km=lai_365(:,:,i);
            albedo_25km=albedo_365(:,:,i);
            wind_25km=wind_365(:,:,i);
            temp_25km=temp_365(:,:,i);
            rh_25km=rh_365(:,:,i);
            pres_25km=pres_365(:,:,i);
            srad_25km=srad_365(:,:,i);
            lrad_25km=lrad_365(:,:,i);
            rn_25km=rn_365(:,:,i);
            sm1_25km=sm1_365(:,:,i);
            lc_25km=lc_365(:,:,i);
            
            c=ndvi_25km>0.011& albedo_25km>0.011& lc_25km>0 & lc_25km<16;
            if numel(albedo_25km(c))>20
                
                albedo_25km_in=albedo_25km(c);
                albedo_25km_in=albedo_25km_in';
                [albedo_25km_in]=mapminmax('apply',albedo_25km_in,PALBEDO);%mapminmax(albedo_25km);
                
                ndvi_25km_in=ndvi_25km(c);
                ndvi_25km_in=ndvi_25km_in';
                [ndvi_25km_in]=mapminmax('apply',ndvi_25km_in,PNDVI);%mapminmax(ndvi_25km);
                
                lai_25km_in=lai_25km(c);
                lai_25km_in=lai_25km_in';
                [lai_25km_in]=mapminmax('apply',lai_25km_in,PNDVI);%mapminmax(ndvi_25km);
                
                wind_25km_in=wind_25km(c);
                wind_25km_in=wind_25km_in';
                [wind_25km_in]=mapminmax('apply',wind_25km_in,Pwind);%mapminmax(wind_25km);
                
                temp_25km_in=temp_25km(c);
                temp_25km_in=temp_25km_in';
                [temp_25km_in]=mapminmax('apply',temp_25km_in,Ptemp);%mapminmax(temp_25km);
                
                rh_25km_in=rh_25km(c);
                rh_25km_in=rh_25km_in';
                [rh_25km_in]=mapminmax('apply',rh_25km_in,Prh);%mapminmax(rh_25km);
                
                pres_25km_in=pres_25km(c);
                pres_25km_in=pres_25km_in';
                [pres_25km_in]=mapminmax('apply',pres_25km_in,Ppres);%mapminmax(pres_25km);
                
                srad_25km_in=srad_25km(c);
                srad_25km_in=srad_25km_in';
                [srad_25km_in]=mapminmax('apply',srad_25km_in,Psrad);%mapminmax(srad_25km);
                
                lrad_25km_in=lrad_25km(c);
                lrad_25km_in=lrad_25km_in';
                [lrad_25km_in]=mapminmax('apply',lrad_25km_in,Plrad);%mapminmax(lrad_25km);
                
                
                sm1_25km_in=sm1_25km(c);
                sm1_25km_in=sm1_25km_in';
                [sm1_25km_in]=mapminmax('apply',sm1_25km_in,Psm1);%mapminmax(sm1_25km);
                
                lc_25km_in=lc_25km(c);
                lc_25km_in=lc_25km_in';
                [lc_25km_in]=mapminmax('apply',lc_25km_in,Plc);%mapminmax(sm1_25km);
                
                dem_25km_in=dem_25km(c);
                dem_25km_in=dem_25km_in';
                [dem_25km_in]=mapminmax('apply',dem_25km_in,Pdem);%mapminmax(dem_25km);
                
                P_25km=[ndvi_25km_in(:) albedo_25km_in(:) wind_25km_in(:) temp_25km_in(:) rh_25km_in(:) ...
                    pres_25km_in(:) srad_25km_in(:) lrad_25km_in(:) dem_25km_in(:) sm1_25km_in(:) lc_25km_in(:)];
                              
                P_25km=P_25km';
                SM_25km_result0=predict(net,P_25km);
                SM_25km_result=SM_25km_result0(1,:);
                SM_25km_result2=SM_25km_result0(2,:);
                SM_25km_result=mapminmax('reverse',SM_25km_result,Pe);
                SM_25km_e(c)=(SM_25km_result>0 & SM_25km_result2>0).*SM_25km_result;
                SM_25km_e=reshape(SM_25km_e,550,600);
                
                SM_25km_result2=SM_25km_result0(2,:);
                SM_25km_result2=mapminmax('reverse',SM_25km_result2,Pt);
                SM_25km_t(c)=(SM_25km_result>0 & SM_25km_result2>0).*SM_25km_result2;
                SM_25km_t=reshape(SM_25km_t,550,600);
                
                SM_25km_et=SM_25km_e+SM_25km_t;%  reshape(SM_25km_et,550,600);
            end
            if i<10
                write_tiff(SM_25km_e,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_e\' ,num2str(i_year),'\', 'E_00' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            elseif i<100
                write_tiff(SM_25km_e,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_e\' ,num2str(i_year),'\', 'E_0' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            else
                write_tiff(SM_25km_e,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_e\' ,num2str(i_year),'\','E_' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            end
            
            if i<10
                write_tiff(SM_25km_t,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_t\' ,num2str(i_year),'\', 'T_00' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            elseif i<100
                write_tiff(SM_25km_t,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_t\' ,num2str(i_year),'\', 'T_0' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            else
                write_tiff(SM_25km_t,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_t\' ,num2str(i_year),'\','T_' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            end
            
            if i<10
                write_tiff(SM_25km_et,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_et\' ,num2str(i_year),'\', 'ET_00' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            elseif i<100
                write_tiff(SM_25km_et,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_et\' ,num2str(i_year),'\', 'ET_0' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            else
                write_tiff(SM_25km_et,['G:\CAS\data\ICTP\ET_Heihe\ET_DNN_TSEB\result_et\' ,num2str(i_year),'\','ET_' num2str(i) 'dds' num2str(j_ns) '_' num2str(j_nl) '.tif']);
            end
        end
    end
end

%%  xlswrite(['G:\CAS\data\ICTP\ET_Heihe\res\a\',num2str(i_year),'\',RR.xls'], R_result, 'Sheet1');


function write_tiff(data,inputpath)
t = Tiff(inputpath,'w');
% Ӱ���С��Ϣ��������Ƚϼ򵥣�
tagstruct.ImageLength = size(data,1); % Ӱ��ĳ���
tagstruct.ImageWidth = size(data,2);  % Ӱ��Ŀ���

% ��ɫ�ռ���ͷ�ʽ����ϸ������3.1��
tagstruct.Photometric = 1;

% ÿ�����ص���ֵλ����singleΪ�����ȸ����ͣ�����32ΪϵͳΪ32
tagstruct.BitsPerSample = 64;
% ÿ�����صĲ��θ�����һ��ͼ��Ϊ1��3�����Ƕ���ң��Ӱ����ڶ���������Գ�������3
tagstruct.SamplesPerPixel = 1;
tagstruct.RowsPerStrip = 16;
tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
% ��ʾ����Ӱ�������
tagstruct.Software = 'MATLAB';
% ��ʾ���������͵Ľ���
tagstruct.SampleFormat = 3;
% ����Tiff�����tag
t.setTag(tagstruct);

% ��׼����ͷ�ļ�����ʼд����
t.write(data);
% �ر�Ӱ��
t.close;
end